





DETERMINATION OF TURBULENT FLOW FIELD OVEI 
A FLAT PLATE USING ENSEMBLE-AVERAGED 
NAVIER-STOKES EQUATIONS 


A thesis submitted 

in partial fulfilment of the requirements 
for the degree of 

Master of Technology 


by 

S.P.Chandra Sekhar 


to the 


DEPARTMENT OF AEROSPACE ENGINEERING 

INDIAN INSTITUTE OF TECHNOLOGY-KANPUR 

April-1994 



AE- Wh- ev- qe ic- per 

2 5 APR ^4 
■ 

' A !M 

lit. ifo. A. I 


CERTIFICATE 


It is certified that the work contained in the thesis entitled 
"DETERMINATION OF TURBULENT FLOW FIELD OVER A FLAT 
PLATE USING ENSEMBLE-AVERAGED NAVIER-STOKES EQUATIONS" 
by S.P.Chandra Sekhar, has been carried out under my supervision and 
that this work has not been submitted else where for a degree. 


O V ^ v- av ' v'-. o. 


Dr. O.P.Sharma 
Professor 

Dept, of Aerospace Engg. 
Indian Institute of Technology 
Kanpur 


Date : 6-4-1994 
Place :I.I.T., Kanpur 



Dedicated to 


my parents 

Sri. S.P. Veera Bhadra Sastry 
Smt. Swarajya Lakshmi 

and to 
my teacher 
Dr. O.P.Sharma 



Acknowledgement. 


First of all I wish to express my deep gratitude to Dr. O.P.Sharma for his 
valuable guidance and encouragement towards the successful completion of 
this work. 1 really enjoyed a great pleasure while working under him. 

I am thankful to my friends V.Durga. Prasad, N.Srinivas, U.V.Sarma 
G.Srinivas, M.Jaya Krishna., S.Arunkumar, Prashant Deb and Shamsi for their 
immemorable help at all stages of this work. 

I am also thankful to my family members for their encouragement and 
moral support during my graduation. 

At last, but not least, I wish to express my sincere thanks to the Computer 
Center and Central Library people for extending their cooperation. 


Date : 6-4-1994 
Place : I.I.T.Kanpur. 


S.P.C.S. 



Abstract 


The compressible, turbulent flow over a thin flat plate is analysed by solving ensemble- 
averaged Navier-Stokes equations. A computer code is developed based on time-dependent, 
multi-dimensional implicit finite difference scheme to solve the governing partial differen- 
tial equations. The method due to Briley and McDonald involving time linearization 
and alternating-directional implicit technique is used. The discretization of the governing 
equations and the corresponding boundary conditions lead to a system of block tridiago- 
nal form of algebraic equations containing the unknown flow variables at each grid point. 
.Assuming uniform initial conditions and solving these algebraic equations result the solu- 
tion of the flow field at the unknown time level. The process is repeated by marching in 
time direction until the solution attains steady state. A zero-equation Cebeci and Smith 
turbulence model is used to describe the turbulence of the flow. A second order arti- 
ficial dissipation model is used at implicit time level to suppress the spatial oscillations 
developing due to central differencing of the convective terms and a fourth order explicit 
dissipation is used to smoothen the solution at each time step. 

The flow field over the flat plate is analysed starting from the leading edge to 
the turbulent region. In the laminar region, the calculated velocity profile is compared 
with the known Blasius solution and a good agreement is found. Next, in the turbulent 
region, the calculated velocity is compared with Simpson’s experimental data. For various 
input parameters, a comparison is made between the code prediction and an empirical 
correlation between the coefficient of local skin friction and local Reynolds no. along the 
plate length and a reasonably good agreement is found. Thus the general validity of the 
numerical scheme and computer code has been established. 
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Chapter 1 


Introduction 


1.1 General 


Most of the flows encountered in nature and engineering systems are turbulent. These are 
highly unsteady with random fluctuations in both space and time. It presents some of the 
most difficult problems both in fundamental understanding of its behaviour as well as in 
obtaining solutions of realistic situations, many of which are still unresolved. An important 
characteristic of turbulence is the richness of scales associated with eddy motion present 
in the flow. The flow over an aircraft wing, cascade flows, rocket motor internal flows are 
some examples of real life turbulent flows. 

In fluids of low viscosity, the steady laminar flow tends to become unstable at 
high Reynolds numbers. This instability due to small disturbances is an initial step in the 
process whereby a laminar flow undergoes transition and then becomes fully turbulent. 
The description of transition again, is a big difficulty due to lack of understanding of the 
underlying phenomena. 

' The flow field can be analysed either by conducting an experiment or calculating 
theoretically. These two methods have their own advantages and disadvantages. Though 
the experiments will give the most realistic answers for many flow situations, however, the 
cost, sensitivity, accuracy etc. of the equipment and measmdng instruments are the limiting 
factors. Moreover, in performing experiments it is not always possible to simulate the true 
operating conditions of the prototype. Furthermore, there is also a need of refinement of 
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the existing measuring techniques for real complex flows. 

On the other hand, in theoretical approach, the flow field physics is translated 
into mathematical equations resulting, in general, in a set of differential equations. These 
generally reflect conservation of mass, momentum, and energy of the fluid, which when 
solved with appropriate boundary conditions give the solution of the flow field. Unlike 
the laboratory experiments, this approach does not entail any problems such as probe size 
limitations, etc. and information can be obtained about quantities which are difficult or 
impossible to measure in the laboratory such as the instantaneous pressure in the interior 
of the flow, hot flow fields etc.. The general solution can be obtained in some simple cases 
by analytical methods, however, in many cases of practical interest, analytical solutions 
are not possible due to the complex forms of the governing differential equations. These 
types of problems can be handled conveniently to some extent by numerical methods using 
a high speed digital computer. 

The development of equations governing the fluid motion started with Navier in 
1822 and the subsequent contributions by St.Venant, Poisson and Stokes led to the present 
day well known Navier-Stokes equations. These are a set of nonlinear mixed type coupled 
partial differential equations, and are considered to be the corner stones of the Continuum 
Newtonian Fluid Dynamics. These conservation equations are valid for the turbulent flows 
as well, however, the values of the flow variables are to be replaced by their instantaneous 
values since it is a fluctuating flow [ 8 ]. 

One approach to deal with the existence of turbulence is to solve the equations 
for a set of boundary conditions and initial values for a large number of instants of time and 
then to compute the mean values ( known as the ensemble-averaged solutions ) . But this is 
a hopeless undertaking even for a most restricted problem such as the case of incompressible 
flow, because it requires large amount of computation as well as computer memory, which is 
beyond the capacity of the present day computers. In view of these limitations, in another 
way, the solution of the turbulent flow makes use of statistical averaging procedures to 
average the governing equations and then solving these mean flow equations rather than 
averaging over the solutions. The averaging can be either time averaging, mass-weighted 
averaging or ensemble-averaging [ 8 ]. Tliis averaging of the equations give rise to new terms 
and require some approximations to close the system of equations. These additional terms 
are related with mean flow variables through turbulence models and can be interpreted as 



apparent stress gradients and heat flux quantities associated with the turbulent fluctuating 
quantities. 

Different models have been developed with different levels of approximations to 
close the set of mean equations. The models can be generally classified into Eddy Viscosity 
models and Reynolds Stress models [ 16 ]. There are various eddy viscosity models of 
varying sophistication. The simplest model is due to Boussinesq ( 1877 ) in which the 
turbulent stresses in the mean momentum equations are expressed as the product of a 
turbulent viscosity and the mean velocity tensor gradient. The turbulent viscosity also 
called eddy viscosity is proportional to the mean velocity. Further development involves 
the well known Prandtl mixing length concept (1925) in which the eddy viscosity is 
expressed in terms of the mean velocity gradient. This tui'bulent viscosity is a property 
of the flow whereas the molecular viscosity is the property of the fluid and this concept of 
eddy viscosity is phenomenological and has no mathematical basis. The eddy viscosity can 
be prescribed algebraically in terms of the quantities such as turbulent kinetic energy and 
dissipation which are obtained by solving the corresponding partial differential equations. 
Depending on the type and number of equations employed for turbulent quantities, these 
are called algebraic or zero equation, one equation and two equation models. Zero-equation 
models are based on mixing length and eddy viscosity concepts and are being widely used 
in practical engineering applications. The turbulent kinetic energy pde, associated with the 
so-called velocity scale , in case of one-equation model and in case of two-equation mode] 
along with the turbulent kinetic energy pde, the turbulent dissipation pde, associated 
with the so-called length scale , are added to the governing equations through the meat 
variables to determine the turbulent viscosity. The one and two-equation models an 
employed when additional details on turbulence quantities are needed [ 16 ] . 

The other models are Reynolds stress model, algebraic Reynolds stress mode 
and large eddy simulation. Reynolds stress model employs several pde’s for component! 
of the turbulence stress tensor. This is one of the most complex models in the curren 
literature. The Reynolds stress model is under extensive development and is employed ii 
complex flow situations, three-dimensional flows, flows with curvature, rotation, blowinj 
or suction etc. In algebraic Reynolds stress models, approximations are made to reduc 
the governing pde’s of the Reynolds stresses to algebraic equations. These equations ar 
coupled to transport equations of the two-equation models to find the turbulence and floi 
field. In large eddy simulation the time dependent, three-dimensional, large eddy structur 



4 


is resolved . through a numerical solution of the time dependent Navier-S fcokes equations 
and a low level model for the small scale turbulence is employed. This kind of simulation 
is prohibitively expensive and it’s use is presently limited to very simple flows. With slight 
variations in different empirical constants of these various models are reported, suitable to 
specific types of problems. But there is no universal model applicable to all types of flows 
so far . 


To reduce the complexity of the problem, different levels of approximations are 
made to the original Navier-Stokes equations such as boundary layer equations, parabolized 
Navier-Stokes equations, etc. [ 1 ]. In parabolized equations, also called thin layer Navier- 
Stokes equations, the mean momentum equation in normal direction is considered but the 
stream wise momentum diffusion is neglected. In boundary layer equations the normal 
direction mean momentum equations are also neglected. These kinds of simplifications are 
found to be reasonable only in certain situations. The need to solve the full equations is 
still remained for a variety of real life situations like flow at a corner, the flow over a launch 
vehicle, etc.. 

The unsteady compressible Navier-Stokes equations are, in general, a mixed set 
of hyperbolic and parabolic typo equations whereas the unsteady incompressible Navier- 
Stokes equations are a mixed set of elliptic and parabolic type equations. If the unsteady 
terms are dropped from compressible equations, the resulting equations become a mixed set 
of hyperbolic and elliptic type equations. These are difficult, to solve because of the different, 
numerical techniques required for each of them. As a consequence, nearly all successful 
solutions of the compressible Navier-Stokes equations have employed the unsteady form of 
equations. The steady state of solution is obtained by marching the solution in time until 
prescribed accuracy in convergence is achieved. Explicit, implicit and weighted average 
of explicit and implicit techniques are the basic methods to solve these time dependent 
equations. 

In explicit methods, the spatial gradients are evaluated using known conditions 
at the old time level. These methods are simpler, easy to code, and vectorizable. The 
implementation of the boundary conditions is also easy. But the major disadvantage is the 
requirement to satisfy the stability condition as dictated by the CFL criterion and thereby 
large number of time marching steps are needed. In the implicit methods, the spatial 
derivatives are evaluated at the unknown time level. This permits the use of higher time 



steps, so faster convergence can be achieved. But for each time step, solution is obtained 
by solving a set of algebraic equations which needs again a big computation. In the implicit 
methods, due to the large time steps, there will be more truncation error, and results 
relatively inaccurate transient path. So this method is a most attractive if the steady 
state solution is only important rather than the transient solutions. The weighted average 
of explicit and implicit method uses the earlier two concepts with some weightages. For 
example, the Crank-Nicolson method, in which equal weightage is used. 

Direct solutions in these methods involve inversion of a block pentadiagonal 
matrix for 2D problems, and block septadiagonal matrix for 3D problems at each time 
level, which consume large CPU time, so approximate factorization is often employed to 
convert this into a simple block tridiagonal system. Tliis in turn introduces a factorization 
error and may lead to stability problems. So some care should be taken regarding the 
time step while developing the solution. Some of the earlier methods utilized the explicit 
techniques, while later development involved the implicit schemes. 

For solving the fluid dynamics equations, there are various numerical methods 
as summarized below: 

1. Finite Difference Method ( FDM ) : In this method a mesh is superimposed over 
the flow field, with the mesh lines along the coordinate directions. The derivatives 
in the pdes are replaced by their finite difference approximations and the unknown 
variables will be evaluated at each grid point by solving the resultant system of 
algebraic equations [ 1 ]. 

2. Finite Volume Method ( FVM ) : In this method the computational domain is di- 
vided into a network of finite volumes. The equations of fluid dynamics are used 
in integral form. The unknown vector and fluxes are staggered, i.e., the unknown 
vector is required at centroids of finite volumes, while the fluxes are required at the 
surfaces of these volumes. The network of the volumes can be arbitrary and numer- 
ical quadrature is used to replace flux integrals in terms of values of the unknowns 
at the centroids [ 1 ]. 

3. Finite Element Method ( FEM ) : The domain is divided into several finite elements 
which can be triangles, tetrahedrons quadrilaterals or hexaliedral elements etc. A 
piece-wise linear approximation to the unknowns in terms of their values at the nodes 
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of the finite elements is assumed. Minimization of the error gives the solution at the 
nodal [ 14 ]. 

4. Particle In Cell Method ( PIC ) : A mesh is required in this method also. The 
transport terms are simulated by movement of particles from cell to coll. The pressure 
gradient and viscous terms are treated by finite differencing. The PIC method is 
thus a combination of Lagrangian and Eulerian description of fluid flow. Further 
development of this model is Statistical PIC Method. In this method stochastic 
model equivalent to the PDE is used for numerical simulation. This approach is also 
called Monte-Carlo Method [ 19 ]. 

5. Spectral Methods ( SM ) : The unknown is expressed in terms of basis functions 
like Chebychev polynomials or trigonometric functions over a mesh. The amplitudes 
of various modes can be obtained from the data at the mesh points. The space 
derivatives can be then obtained by differentiating the Fourier expansion. Then the 
PDE is reduced to a set of ODE by taking its inner product with basis functions. 
The resultant ODE can be solved by differencing the time derivatives [ 7 ]. 

6. Cellular Automata or Lattice Gas Dynamics ( CA ) : The computational domain 
is divided into a fine mesh called lattice. At the lattice points particles are located 
having velocities with directions along the mesh lines passing through the lattice 
points. The transport term is simulated by shifting the particles to neighbouring 
lattice points. A collision cycle of calculations is then performed dining which pre- 
collisions velocities of particles are replaced by post collision velocities [ 10 ]• 

Out of these methods FDM, FVM, FEM are well developed while the remaining 
methods are presently in developing stage. For a simple geometry, with uniform grids, 
Finite Difference Methods are suitable. These are simple, easy to understand, having 
relatively less algorithm complexity, and recommended by many researchers also. 

Several computer codes have been reported in the literature to solve different 
levels of approximated equations during the past decade or so. Some of them are MINT 
developed by McDonald and his associates, ARC-2D and ARC-3D developed by Pulliam 
and Steger, TEACH developed by Spalding group, NISA/3D Fluid developed by EMRC. 
Troy, USA etc. Few of them are made available on commercial basis also. Relevant journal 
articles appear in AIAA Journal, Computers and Fluids, Computer Methods in Applied 
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Mechanics and Engineering, Journal of Computational Physics, Journal of Fluid Mechan- 
ics, Numerical Heat Transfer, Transactions of ASME., Journal of Fluids Engineering, and 
some of tlxe important books in this area are given in the list of references. 


1.2 Present Investigation 


As mentioned earlier there are different forms of approximated Navier-Stokes equations 
which are applicable in particular situations only. In the cases of mixing of the viscid and 
inviscid flows, chemically reacting flows or flows having complicated geometries such as 
flow at slots, submerged nozzles of rocket motor etc, the approximated equations will not 
represent the entire physics of the problem and in such situations the simplified equations 
are not valid. Therefore, there exists a need to solve the full equations for the entire 
domain. 

In the present investigation, for simplicity, the flow over a thin flat, plate is 
considered. Adjacent to the plate there will be a boundary layer and it will be laminar and 
then gradually undergoes transition and becomes fully turbulent as the flow moves from 
the leading edge to trailing edge of the plate. The coefficient of wall friction, the thickness 
of the boundary layer and other flow variables vary along the plate length. Getting a 
detailed steady solution for a turbulent flow field over a thin flat plate by solving the 
ensemble-averaged Navier-Stokes equations is the goal of the present investigation. 

The steady state solution is obtained by marching in time direction until the 
convergence in time is reached . A mesh, clustered , near the wall is superimposed over the 
flow field and an implicit finite difference method is used to solve the governing conservation 
equations. Briley and McDonald [ 5 ] developed a split linearised block implicit numerical 
scheme to solve the multi-dimensional unsteady Navier-Stokes equations using alternating 
direction implicit technique and tested it for a laminar flow in the entrance region of a 
duct. Further, Sabnis et al. [ 21 ] used it to solve the solid propellent rocket motor internal 
flows. The same procedure is adopted here to solve the external turbulent flow over a flat; 
plate. This method is a one-step method, as opposed to a predictor-corrector method, 
and requires no iterations to compute the solution at single time step. Cebeci and S mi th; 
turbulence model [ 8 ] is used for the turbulent viscosity in turbulent boundary layer. 
Further a second order artificial dissipation model and a fourth order smoot hin g terms are 
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included to suppress the spatial wiggles. A computer code is developed to solve this flow 
field and validated with the available data. 

The details of the method of solution , the results and conclusions are described 
in the following chapters. 



Chapter 2 


Problem Formulation 


2.1 The Flow Field Description 

The flow over a thin flat plate is analysed by solving the general Navier-Stokes equations 
for a turbulent, flow. A schematic diagram of the flow field is shown in Figure 2.1. It is 
considered that the fluid is a viscous, compressible, chemically non-reacting, perfect gas 
with zero bulk viscosity coefficient, constant molecular viscosity, thermal conductivity and 
specific heats. The flow over the plate is considered with zero pressure gradient. The 
plate is very thin with smooth surface and is fixed so that the flow field oscillations due to 
turbulence are not transmitted to it. Thermodynamically, the plate is assumed as perfectly 
adiabatic and no heat transfer takes place. Alternate conditions for wall thermal conditions 
and for pressure gradient can also be incorporated. 

Adjacent to the plate, there will be viscous and thermal boundary layers. Near 
the leading edge, the flow in the boundary layer is laminar and it becomes unstable as 
it moves in downstream. Then the flow undergoes transition in which the disturbances 
grow, and these disturbances increase in the flow and becomes fully turbulent as it moves 
further downstream ( see Figure 2.1 ). The flow in the turbulent layer is more dissipative 
in nature, and no longer be similar. The boundary layer thickness, the velocity gradients 
at wall surface, etc. vary along the length of the plate. There will be heating of the fluid 
also due to the viscous effects as it moves downstream. 

The fluid motion at any instant of time is described by the well-known Navier- 
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Stokes equations, which are the statements of conservation of mass, momentiun, and energy 
of the fluid. In order to deal with turbulent How, as mentioned earlier, the mean flow 
equations can be obtained by using ensemble-averaging. The ensemble-average of a variable 
is defined as the arithmetic mean measured over the many identical experiments [ 17 ]. 

Denoting the ensemble average of a fluctuating variable u(t) as u(t), it is defined 
as 

u (*) = jj E [ \i (2- 1 ) 

Similarly, the ensemble average of a function g(u), say, g(u) is defined as 

g(u) = jijm [«(*)] Ij ( 2 - 2 ) 

~ >0 ° v j= 1 

An instantaneous variable, say, u can be represented as a sum of a mean quantity u and a 
fluctuating quantity u' . Thus 

u = u + u (2.3) 

In the process of averaging the following rules have been used for averaging different 
quantities. 


uu' 

= U 'll = 

0 

(2.4) 

‘Hi Uj 

= ( Ui + 1 

J,'i) (uj -1- Uj ) 



= Ui Ui + 

■u'i u'j 

(2.5) 


and so on. The application of this averaging procedure to the governing differential equa- 
tions give rise to the mean equations. The resulting mean equations, thus contain some 
additional terms in view of the turbulent fluctuations and are modelled by introducing the 
concept of turbulent viscosity and conductivity as described in the previous chapter. 


2.2 Governing Equations 

The governing equations of a turbulent flow field over a thin adiabatic flat plate, such 
as the ensemble-averaged Navier-Stokes equations, describing the conservation of mass, 
momentiun, and energy of the fluid along with the equation of state for perfect gas are 
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given below in vector form [ 22 ]. The variables in the following equations are ensemble 
averaged physical quantities, and for simplicity the symbols overbar and tilde are removed. 

Continuity equation : 

fj + V-Utf) = 0 (2.6) 

Momentum Conservation equation : 

j t (pU ) + V-(pUU) = -Vp + V-r (2.7) 

Where 


r is the stress tensor ( both molecular and turbulent ) given by 

Tij = 2 [i e jj &ij -| n e ffV-U6ij 

eij is the rate of strain tensor, given by 

1 / dUj dUj \ 

C,; 2 \ dxj dxi ) 

H e jj is the effective viscosity and is the sum of the molecular and turbulent 
viscosities. That is, 


( 2 . 8 ) 


(2.9) 


»eff = A* + to (2-10) 

The turbulent viscosity is obtained from the turbulence model and will be described later. 

Energy conservation equation : 

f t (ph) + V-(pUh) = + * (2.11) 

Where 


$ is the viscous dissipation per unit volume, is given by 


$ = 


'*<# 


2^ey - - (V-Uf 


( 2 . 12 ) 



q is the heat flux vector, is given by 
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q = - k c ff^ T (2.13) 

k c jj is the effective conductivity and is the sum of the molecular and turbulent 
conductivities. That is, 

k e jj = k + h (2-14) 

The turbulent conductivity is obtained from the turbulence model. 

Equation of state for perfect gas : 

p = pRT (2.15) 

The pressure term in the above equations can be eliminated by using the equation 
of state for perfect gas. The energy equation Eq.(2.11), can be rewritten in terms of 
temperature as follows. 

C v ~{pT) + V • ( p t7 T ) = - V -q + $ - pRTV -U (2.16) 

These governing equations are coupled, nonlinear and mixed type partial differ- 
ential equations. The momentum equation written in vector form Eq.(2.7) leads to two 
momentum equations in each x and y directions- There are four pdes such as Eqs.(2.5), 
(2.7) and (2.16) and an equation of state for perfect gas Eq.(2.15) and five unknown vari- 
ables such as p, U ( consisting of u and v ), T and p and these can be solved with proper 
boundary conditions as described below. 

2.3 Boundary Conditions 

The following boundary conditions can be obtained by studying the flow field over the flat 
plate. The above set of equations are solved subjected to these boundary conditions. 

Along the plate surface no slip conditions hold good and the plate is assumed as 
adiabatic, thus 


1 


0 


= 0 


(2.17 

(2.18; 


Along the wall: at y = 0 : 

u — v 
OT 
dy 

The density at wall is calculated from the simplified form of the continuity equation al 
wall surface. 

si n s) 

(2.19) 


dp 0 

m + = 0 


Pressure gradient along the dominant flow direction is taken as zero. This gives 
the boundary conditions at the free stream boundary that the density, velocity, tempera- 
ture are constant and equal to free stream values. As this boundary is taken at sufficiently 
far from the plate surface, the gradients of the density, velocity, temperature in normal to 


plate at this boundary are set to 

zero. 

Thus 



Along the free stream boundary 

: as y 

— 

oo : 


u = ; 

v = 

0; 

T = 

-- ; 

dp _ 

da 


dv 

d T 

dy 

dy 


dy 

dy 


Poe 

0 


( 2 . 20 ) 

( 2 . 21 ) 


At the upstream boundary the flow is entering uniformly at zero angle of attack. 
This condition leads to 


Along the upstream boundary : at x — 0 : 

u = We*, ; u = 0 ; T = Too : p = Poo (2.22) 


At the downstream end of the plate, there is lack of information regarding the 
flow, it is, therefore, necessary to model the boundary conditions at this end. This is done 
by setting the second derivatives of the velocity components, temperature in x direction 
are zero. The boundary conditions for density is modelled through equation of state for 
perfect gas. Thus 


Along the down stream boundary : at x = L : 

d 2 u d 2 v d 2 T 

dx 2 dx 2 dx 2 


(2.23) 


For the initial conditions, uniform flow is assumed in the entire domain, and the 
velocity components at the wall are set to zero. 
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2.4 Non-Dimensionalization 

The governing fluid dynamic equations are very often put in the non-dimensional form as 
it leads to some advantages. For instance, that the flow variables can be normalised so 
that their values will fall in between certain prescribed limits such as 0 and 1. Also the 
characteristic parameters of the flow field can be obtained such as Mach number ( M ), 
Reynolds number ( Re ), Prandtl number ( Pr ) which can be varied independently. A 
generalized way of non-dimensionalization is used as described below. 

For non-dimensionalization, the various characteristic parameters are taken as 
follows. The plate length ( L ) is taken as characteristic length ( L r ), the free stream den- 
sity, velocity, temperature, molecular viscosity, molecular conductivity ( u^, X^, /i x . 

are taken as corresponding characteristic quantities ( p r , u r , T r , p. r , k r ) and for time, 
the convective time ( L^/u^ ) is taken as characteristic time ( t r ). And the non- 
dimensionalization procedure is as follows. 


x y 

V y = V 

P = 

P_ 

Pr 

(2.24) 

U V 


f 


U = — ; V = — ; 

u r u r 

T = 

T r 

(2.25) 

_ p . f, ~ h . 

t = 

i 

(2.26) 

PrU; ’ p r u; ’ 

Lrju r 

A , 

p = — -; k = 

Pr 

k 

K 


(2.27) 


Corresponding to these non-dimensional expressions, different first and second order deriva- 
tives are expressed and substituting these in the governing equations Eqs.(2.6 to 2.16), 
various non-dimensional parameters are obtained such as : 


Mach number 

= M = 

Ui r 

C 

(2.28) 

Reynolds number 

= Re = 

Pr * Ly 

Pr 

(2.29) 

Prandtl number 

- Pr = 

P r Cp 
l. 

(2.30) 


k ’7‘ 


After non-dimensionalization the governing conservation equations are rewritten in a 
natrix form as follows. 


OHjp) 

dt 


D (b) + 5(b) 


(2.31) 
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Where <f> is a column vector, H and S are column vector functions of <p and D is 
a column vector whose elements are spatial differential operators. The various elements of 
these matrices are given below. 


W) 


( p, pu, pv, pT ) 


(2.32) 


DM) 


( P u ) 


JL 

dx 

ai ( pu + Ttx-' 


d , , , Veff d 2 v 

ei (puv) + & a? 


P eff d 2 u 4 du 

/2e cbr 2 3 Re dx dx 

2 off du 
3 <9y dx 


£(PUT) + (1 - 7)pT^ 


+ 


lk 


iff 


d 2 T 


+ 


7 OT Ok, 


if 


Re Pr dx 2 Re Pr dx dx 


(2.33) 


Dy{4>) 


/ 


( ) 


d_ 
dy 

JL ( mv) + _ j_^dv 

dy V J Re dy 2 3 i?.e 5.x <9;y 


5 / 2 1 m 

<%( pv + 7M 2pT 


P eff d 2 v ^ e jf ^ 

i?.e dy 2 3 i?.e dy dy 


d . . 

^(p C r ) + (i- 7) ,r- 


' TT> Tf^ rv I 


7 <yT dk 


iff 


RePr dy 2 Re Pr dy dy 


(2.34) 


J 



1 « 




3 Re dx ■ J 0 y l c hj 


iMJL (v . a)+ Tsff(^ + 

3 Re dy K 1 dx \ 6 hj 
7(7 - 1) M 2 ^ [*] 


fltA 

dx) 

dv\ 
dx j 


J 


2.5 Non-Dimensional Boundary Conditions 


(2.35) 


The following non-dimensional boundary conditions are obtained by normalising the bound- 
ary conditions given in Section 2.3. 


At y = 0 : 


u = v = 0 
DT 


dy 


= 0 


(2.30) 

(2.37) 


The density at wall is calculated by simplifying the continuity equation at wall as follows 

= 0 (2.38) 


dp d 

at + ¥ y {pv) 


As y — > oo : 


At x = 0 


u = 

i; 

v = 0 ; 

T = 

i; 

dp 


du 

dv 

dr 

dy 


dy 

dy 

dy 

u = 

1; 

v = 0 ; 

T = 

i; 


At x 


L : 


dru 

dx 2 


d 2 v 

dx 2 


d 2 T 
dx 2 


P = 1 

= 0 

P = 1 

0 


(2.39) 

(2.40) 

(2.41) 

(2.42) 


Similarly, the initial conditions in the entire domain become unity, and the ve- 
locity components at the wall are become zero. 
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2.6 Turbulence Modelling 

In order to close the system of mean flow equations, the turbulent stresses and turbulent, 
heat flux vectors have to be modelled through the mean flow variables. Therefore, following 
the laminar flow description, the turbulent stress is expressed as a product of a turbulent, 
viscosity and mean flow velocity gradient and the turbulent heat flux as a product of 
turbulent conductivity and temperature gradient. As mentioned earlier different models 
have been developed to describe the turbulent viscosity with different level of empiricism. 
Here the so-called Zero- equation model developed by Cebeci and Smith [ 8 ] is adopted. 

This model is based on the eddy viscosity and mixing length concepts. The 
turbulent boundary layer is regarded as a composite layer of a thin laminar sub-layer, an 
inner layer and an outer layer. The thin laminar sub-layer is neglected in the present 
model. The inner region is much smaller than the outer region, with thickness about 10 
to 20 % of the total boundary layer thickness. In the inner region, turbulent viscosity 
is expressed algebraically in terms of the mean velocity gradient and a mixing length . 
The mixing length is expressed in terms of distance from the wall as earlier proposed by 
VanDriest ( 1956 ) [ 15 ] and subsequently extended by including a damped law of wall. In 
the outer region, turbulent viscosity is described in terms of the velocity defect as proposed 
by Coles [ 8 ]. 

In most practical calculations, it is necessary to deal with the laminar, transient 
and turbulent domains of the boundary layer. According to Emmon’s spot theory ( 1951 ) 

[ 13 ], the flow in the transition region is intermittently turbulent. The transition region 
can be accounted by considering the expression for intermittency factor given by Chen 
and Thyson [ 9 ]. This expression is based on Dhawan and Narasimha’s intermittency 
expression ( 1958 )[ 11 ] for incompressible flows. 

The various expressions of the Zero-equation model due to Cebeci and Smith are 
given below. 

In the inner region : 0 < y < y c : 

o du 

H = P l 2 V Itr 

oy 

= * V [ 1 - exp ( - y/A ) ] 


where l 


(2.43) 

(2.44) 
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In the outer region : y c < y < 6 : 

Pt = ft a 


L < 


u ) dy 


7 tr 


(2.45) 


y c is determined from the continuity of the eddy viscosity across the boundary layer. The 
inner region expression is applied outwards from the wall until the two expressions give the 
same value of the turbulent viscosity ( p, t ) The variables appearing in the above equations 
are as follows. 


A is damping constant given by 

A 

A + = 

K = 

For flows with no mass transfer : 

N = 

For flows with mass transfer : 

\ 2 


A 


+ V 


Tu, 


ft N \ p w 


ft_ 

Pw 


26 

0.41 


1 - 11.8 


t^iv 


- 1 


Pw 


p 


N = 


I Poo 
Poo V Pxu 


p_ 

v:\ 


exp 11.8 v. 


> Pw 


P 


+ exp 11.8 v 


i ,„+ /I w 


P 


where 


V 


P w,x) du 0 


u r 


P 


p ui* Ox 
1/2 


V. 


+ _ 


U r 


For values of R 0 < 5000 : 

a = 

n 

Zi = 


o.om^L 
i + n 

strength of the wake 

0.55 1 - exp ( - 0.243 z\^ 

Rq _ - 
425 


0.298 zi ) ' 


(2.46) 

(2.47) 

(2.48) 


(2.49) 


(2.50) 

(2.51) 

(2.52) 

(2.53) 

(2.54) 


(2.55) 

(2.56) 
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For values of R 0 > 5000 : 


a = 0.0168 (2.57) 


The intermittencey factor, y tr accounts for the transition region. For two dimensional 
flows, it is defined as 


1 

. & 

1 

H 

II 

8" 

G ( x — Xtr ) ( 1 — dx 'j 

\ J Xf r U'qq J 

(2.58) 

Xfr is the x location where the transition begins. G is an empirical factor given by 


G = 

— - SSL f> P 1 34 

C 2 v 2 Xlr 

(2.59) 

R^Xtr ~ 

Poo Woo *tr 

(2.60) 


/^OO 

c = 

60 + 4.86M 1 ' 92 

(2.61) 


For 0 < M < 5 


For a two dimensional flows, the transition Reynolds number can be estimated from the 
following empirical expression given by Cebeci and Smith [ 8 ]. 


Re 0(r 


1.174 


{ 22400 

l 


Re 


, 0.46 

■X lr 


(2.62) 


A constant turbulent Prandtl number equal to 0.92 is assumed ( for air ) and is used to 
calculate the turbulent conductivity. The expression for turbulent Prandtl number is given 
below. 


Pr t = = 0.92 (2.63) 

Here the turbulent kinetic energy is not used in the turbulent viscosity formulation, it 
implies that the mean motion is uneffected by turbulence intensity and the length scale. 


This model can be used in two dimensional flows with mild pressure gradients and 
mild curvature with no flow separation or rotation effects. It is useful in computing three 
dimensional boundary layers with small cross flows and mild pressure gradients having no 
curvature or rotation effects. Lakshminarayana [ 16 ] concluded that for an attached two 
dim ensi onal boundary layer with a mild pressure gradient, there is very little advantage 
in restoring to a hig h er order models. The algebraic model is adequate for the prediction 
of the mean velocity field. However, for flow with strong adverse pressure gradient ( close 
to separation ) and flows with -strong favorable pressure gradient ( laminari/ing ), this 
formulation is only adequate but not sufficient. 
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2.7 Artificial Dissipation 


The steady flow over a flat plate is generally a high Reynolds number flow except in the 
region near the leading edge. Central differencing of the first derivatives in the Navier- 
Stokes equations equations leads to spatial oscillations and further the solution diverges 
also some times at these liigh Reynolds numbers. To suppress these spatial oscillations, 
excess viscosity is provided, while developing the solution, through artificial dissipation 
models. This artificial dissipation can be provided either by adding the artificial viscosity 
terms explicitly to the governing equations or one side discretization of the convective 
terms such as upwinding or hybrid schemes [ 19 ]. Various models of adding artificial 
dissipation for high Reynolds flows have been investigated by Shamroth et al.[ 23 ] and 
evaluated in the context of a one dimensional model problem. 


In the present application, a term of a form 


_ ( 21 \ 

dx y aH dx ) 


(2.G4) 


s added in the unknown time level to the governing equations [ 21 ]. This explicit dissipa- 
ion is added in the dominant flow direction only ( that is in ^-direction ), but not in the 
lormal direction ( y -direction ), so that the boundary layer effects will not be affected. The 
ariable (j> denotes the velocity component u for :c-direction momentum equation, velocity 
omponent v for y-direction momentum equation, density p for continuity equation and 
he enthalpy h for energy equation. 


The coefficient p art is obtained from an expression given by 


p u Ax 

where 

— ( ) [f 1 + A l art } 

' & d ' 

(2.65) 

Ax 

= grid spacing at the point in question 

(2.66) 

7 * 

= 0 for continuity equation 

(2.67) 


= l l rjf f° r momentum equations 

(2.68) 


11 off 

= p for energy equation 

(2.69) 


re parameter cr d initially taken as 0.5, which corresponds to the well-known concept of 
dl-Reynolds number of 2. After obtaining a converged solution with = 0.5, its value 
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is reduced to 0.05 and further to 0.005, and the solution procedure further marches in time 
to obtain the final converged solution. 


In order to suppress the high frequency oscillations in the solution a fourth order 
explicit dissipation term is added at known time level to the right hand side of the governing 
equations as follows. 


- e e 


(Ay)' 4 


dy 4 


(2.70) 


This fourth order term does not effect the formal accuracy of the algorithm. 
The — ve sign to the above fourth order derivative is required in order to produce positive 
damping. The smoothing coefficient e c should be less than approximately 1/16 for stability 
[ 2 ], presently 0.05 is used. 


2.8 Non-Uniform Grid Stretching 

Local refinement of mesh where there will be large gradients is always advisable to predict 
them more accurately. In the present problem, the velocity components, temperature, and 
density vary largely near the wall. So it is better to have a small grid spacing near the wall 
and large grid spacing in the outer region rather than uniform small grid spacing in the 
entire domain. But discretization of the governing equations with variable grid spacing is a 
difficult procedure. This variable grid spacing can be handled conveniently by transforming 
the coordinates ( x, y ) to another set of coordinates ( ?]. ( ) with uniform grid spacing. An 
analytical coordinate transformation devised by Robert’s [ 20 ] is adopted here for a non 
uniform grid generation and is described below. 


If N y grid points are to be used in the range of 0 < y < 1 and if steep gradients 
are to be expected in a region of thickness of a near y = 0, then Roberts transformation 
<(y) is given by 


C( y ) = N, + (N v - l)log 

where 


y + a — 1 


y + (i + 1 
a 2 = 1/(1 — a ) 

Then the corresponding partial derivatives are as follows. 
d d( d 


/log 


a — 1 
(i -f- 1 


dy 


% <yc 


- * h t \ 1 

i h r„ kanruh 

„ & livr/of 


(2.71) 

(2.72) 

(2.73) 





where b 


= ( - + 1 
V a - 1 


*v-l 


(2.77 

The transformation in cc-direction is ontional ... 
transition and turbulent boundary layers and the ’ CC ° rdlng to the lei igths of lamina: 
zoaes ' The same transformation can be uTed ta T“ ° f S ° Iuti ™ » differen 
corresponding to a plane of uniform grid spacing f „ r alS °' The ( *■ » ) Plan 

stretching on N, and <7 are shown in Figure 2.2 and 2 3 ? ^ aDd depeudeQ<: y of grit 



Ymax= 01 


! 

i 

i 

i 

■ 

Free Stream 
Conditions 
( Entering ) 


21 

c 

o 





















ts» 





















II 

|i 16 

o 

z 

ii 

6 








































j 












































i 















i 





1 

























" j 
















: - 

j— - 

— 

— 

j 

— 

— 

— 

— 

— 

— 

— 

— 

— 


— 

— 

— 



i 

. . .... 

t 

zm 

zzz 

i 

— 

— 

— 

— 

— 

— 

— 



— 


— 

— 

— 





P — . 

— . — 



— 

— 

. 

— — 

— 

— ..... 



~~ 


- . ■ 

■ ' 

... .... 

I*'" 

— - — - 

r- — • 




— 3rd 

~-rrZ 

xr:: 



=v ; ii 


| rUjr 

=: ”■ ' * 






r.~ i 



slzlz 







1 i ^ 6 ii 16 2i 

Xndrs= 0.0 Xmax= 10 

Ymin= 0.0 


Flow Direction 

( no of grids ) 


Figure 2.2: Non-Uniform Mesh. For N x = N,, = 21. a x = 0.99, cr y = q.O 1. 
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Figure 2.3: Dependency of Grid Points Cluster on Various Parameters 



Chapter 3 


Method Of Solution 


3.1 General Outline 


The governing equations of a turbulent flow over a thin flat plate and the corresponding 
boundary conditions along with the required turbulence modelling and artificial dissipa- 
tion have been described in the previous chapter. An alternating-directional implicit finite 
difference method is used to solve these set of coupled partial differential equations numer- 
ically. 

In the governing equations, the time derivatives arc replaced by a forward time 
difference approximation. The terms involving the non-linearities at the implicit time 
level are linearized by expanding them about the known time solution by using Taylor’s 
expansion. The spatial derivatives are replaced by the central difference approximations at 
implicit time level, for both first and second order derivatives. An alternating-directional 
implicit scheme is used to reduce the complexity of the computing the solution. This 
converts the two dimensional problem into a set of one dimensional problems corresponding 
to each line of grid points of the mesh. 

This leads to a system of coupled linear difference equations having a tridiagonal 
block-banded structure. By solving these block-banded equations the unknown variables 
at n + 1 time level can be obtained. 

Initially, turbulence viscosity is set equal to zero for a few time marching steps. 
The turbulence viscosity is then calculated from the turbulence model as discussed earlier 
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by using the mean flow solution over the flow field at the known time level. This turbulent 
viscosity is used for solution at the unknown time level. Thus the inarching process is 
repeated until the solution attains the steady state. As the solution approaches asymp- 
totically to the steady state, the marching procedure is stopped when tin 1 consecutive 
solution’s difference falls within the specific minimum value, in the present, analysis, this 
is taken as 0.000001. 

The various steps involved in the solution procedure are discussed in detail in 
the following sections. 

3.2 Time-Linearization 

In order to linearize the non-linear terms appearing in the implicit time governing equa- 
tions, various linearizing techniques are available in the literature. Out of them the method 
of extrapolating the coefficients, is adopted here. This method is having certain advan- 
tages, but it makes the equations complicated. This procedure requires no iterations to 
update the coefficients at any time level so that the computational cost will be reduced 
comparatively. It may be pointed out that rather than utlizing the computation time, 
as in other linearizing methods, to iterate only for the sake of reducing the truncation 
error associated with linearization, the overall accuracy might be improved in the present 
method by using the same computation time to reduce the marching time step size. 

The technique is based on Taylor’s expansion of non-linear implicit terms about, 
the solution at the known time level t n . This leads to a one step, two level scheme, which 
is being linear in unknown quantities and can be solved efficiently with out any iteration. 
Formally the truncation error of this procedure can be made as small as required. Here, 
in view of computational cost, first order approximation is considered, wliich neglects the 
second and higher order terms in Taylor’s expansion [ 1 ] . The procedure of linearization 
is explained below. 

Considering the Eq.(2.31) at implicit time level 

1 


dt 


D(0)" +1 + S (?) 


(3.1) 



Replacing with implicit time difference approximations, 
( H - H \ 


£>(?;)“" + s(d>)" 


At 


The various quantities are expanded as follows. 


71+1 


(3.2) 


H (<£) n+1 
D ( (f>) n+1 
S (4>) n+1 


H^r + 


D (cj>T + 


s(t) n + 


SHWY f^) &t+0(At y 


dt 


dt ) 


d£wy (a# 

dt J \ dt 

dsj^y ( dt' 

d<t> ) \ dt , 


At + O (At) 


At + O (At y 


(3.3) 

(3.4) 

(3.5) 


The term ( ^ ) At is approximated further by ( <j) n+1 — (j )' 1 ) . Thus the Eq.(3.2) be- 
comes a linear implicit time difference form as follows. 

dH(t) n ( t n+1 - 4> n \ 


dt 


At 


= D(t) n + s(ty + 

f OD (t) n dS (<l))' 1 \ ( t n+1 ~ t n 


dt 


0(h 


At 


(3-6) 


On examination, it can be seen that Eq.(3.0) is linear in the quantity ( <f> n+l - t n ). 
Computationally, it is convenient to solve the Eq.(3.6) for ( c/> nM — (j) n ) rather than cj> n+l . 
This also simplifies the roundoff errors. So introducing a new variable $ such that 


$ 

so that 'E ,n+1 

and T” 

Rewriting the Eq.(3.6) in a convenient way 

(A + At B) $ n+1 

where A = 

B 


= t - 

t n 


(3.7) 

= t n+1 

— 

t n 

(3.8) 

= 0 



(3.9) 

At [ D (t) n 

+ 

s (t) n } 

(3,10) 

0H(tT 

dcj, 

At 

( os(t) H ) 

{ dt J 

(3.11) 

dD(t) n 

dcj) 



(3.12) 


Applying the above linearization procedure to each of the conservation equations 
the following forms are obtained. 
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Continuity Equation : 

The continuity equation in implicit time difference form, is 

f p n+1 - p n \ 8 n+ , 8 


tl |- 1 


At 


0 


( 3 . 13 ) 


the nonlinear quantities by Taylor’s expansion as explained before 

At. + O ( At )~ 


Expanding the 
( pa r +i = ( pu r + 

= ( pu ) n + 


&' 

x dp 

m + 

d 

du 

( pu ) 

du 

Of. _ 

n 

At 

“*( 

' p n+1 - p n ' 

) + 

f. 

u n\ 1 

- U n 

At 


At 


= ( pu ) n + u n '£" +1 + p’ 1 Tf 1 


At + O (At) 2 
(3.14) 


Similarly, 
substituting the above expressions in 




' P 

- At 


( pv ) n+1 = ( pv ) n + + p n <2 

Eq.(3.13), the continuity equation becomes 
d ' "’” A1 + At ~{ v n ^; +1 + p n *': +1 ) 


(3.15) 


^" +1 + A t 4- ( U n K +1 + P' l K I . — ' n 

(It. \ * J Qy 


ax'. - 

a . a . 

ad'“ ) + -(pv) 


(3-K5) 


x- momentum equation : 

The x-momentum equation is written in an implicit time difference form. For the sake of 
convenience, the cross derivatives in the viscous dissipation term are treated at explicit 
time level, ( that is at known time level ). This makes the splitting of equations for 
alternating direction is easier. 


( P u ) 


n- 1-1 


( pu y 


d / ,v>+i 8 , 


d 


At 


\n-fl 


+ 


9p n+1 , M 

dx Re 

J_ lAsff 

Re dx 


dx 

d ‘) H 4*1 


dy 




<^9 n-f*l 

o~u 


4 duy i+l 
3 dx 


dx 2 dy 2 

2 dv n+1 ^ 


+ 


1 d 


3 dy 


+ 


3 dx 
d ' l cff ( 


du' 


dy V dy 


du 11 

dx 

+ 


+ 



dv n 

dx 


(3.17) 


The various nonlinear terms are expressed as follows 

(pa 2 ) n+1 = (pu) n + (u 2 ) n V* +1 + (2pu)"t;: +1 (3.18) 

( puu ) n+1 = ( puv) n + (pu) M 'l>" 11 + ( pu )"■'!'" 11 + ( au ) n 'Pp l 1 (3.19) 
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and like wise for other terms. Further, the term pressure p is replaced by pT/'yM ", by 
using the equation of state for perfect gas. The x-momentum equation in the linearised 
implicit form becomes 


+ 

+ 


u n r +1 + 1 


+ At. 


+ 2 (,,«)” n ) 


At 
At 

7 M 2 

At 
Re 

- At 


d 


— ( ( ^ )"*” +i + ( pu rv;, +i + ( pv m +i ) 


( p n w ^ +1 + T n ^" +1 


dy 

JL 

dx 

^eff f 4 9$» +1 _ 2 
9s 3 9s 3 9;y 

d , *\n 1 9pT n 


i'vff 

Re 


+ 
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dx 

L 

At '‘eg 
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( ) + r«75 n. 


Re 
At 


d 9 n rV) 

"U o~u 

dx 2 9y 2 


7 M 2 9s 

+ • 
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dpuv n 
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tV> 1 
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dx 2 
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&P e ff /4 9u’ 1 _ 2 dv‘ 
dx l 3 dx 3 dy 


d 2 ^l +i 

— + 

dx 2 




d 2 ^\ 


dy 2 


d 2 v n 

dxdy 


As f ay ay 

dy l dy dx 


(3.20) 


y-momentum equation : 

Similarly, the y-momentum equation is written in implicit time difference form, again the 
cross derivatives are written at old time level. 


( pv ) n+1 — ( pv ) 


At 

dp n+l 
dy 


+ 


+ 


Re 


As. ( Su 

dx l dy 



(3.21) 


Writing in the linearised form, after expanding the non-linear terms according to Taylor’s 
expansion as described before and replacing the pressure term p by pT/yM 2 
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1 + (i>»rK n 
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At ii eff 
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(3.22) 


Energy equation : 

The energy equation in the implicit time difference form is given below. To reduce the 
complexity, the dissipation term $ is treated explicitly, though it can be discretised at 
implicit level. 


(1) 
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( pt y 
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n-f 1 


+ 
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dy 
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dT 9 jAff dT Q^eff 

dx dx dy dy 


(3.23) 


Where $ is the dissipation function, defined earlier in Eq.(2.12). The pressure term p is 
replaced by pT/yM 2 , as before. Writing the energy equation in the linearised form, after 
expanding the non-linear terms according to Taylor’s expansion as described earlier, 
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3.3 Alternating-Directional Implicit Technique 


Direct solving of the all coupled equations in the entire domain needs inversion of a penta- 
diagonal block banded matrix ( for 2-D problems ). This needs large memory as well as 
large computation. These can be reduced by applying an alternating-direction implicit 
technique, which results an easily solvable narrow ( tridiagonal ) block banded matrix. An 
alternating-direction implicit technique, developed by Douglas and Gunn [ 12 ] is adopted 
here. 

The operator D in Eq.(2.31) is assumed that it contains only the derivatives of 
first and second order with respect to# and y, but no mixed derivatives. This D operator 
can be split into two operators, D x and D y associated with the x and y coordinate directions 
and each having a functional form 

D ‘ = ( - 3 ' 25 ^ 

* - F (*-ww) (3 ' 2(i) 

Then the Eq.(3.10) becomes 

[A + At ( B x + B y )} ’3/ n+1 = A t \{D X + D y )r + S n ] (3.27) 

The Douglus and Gunn representation of the above equations can be written as a two step 
solution procedure as follows. 

First ADI step : 

[A + At (J3,)] V* = At [ (D x + D y ) (IP + S n ] (3.28) 


Second ADI step : 


[A + At (By) ] 4/ nM 


A T 


(3.29) 
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where '£* is an intermediate solution. 


The major attraction of this method is that the intermediate solution is a consis- 
tent approximation. Furthermore, the steady solution when T' 1 = vj>* = T' H1 = 0 

; is independent of At . Thus, physical boundary conditions for T' l+I can be used in the 
intermediate steps without a serious loss in accuracy, with no loss for steady solution [ 5 ]. 
In this respect the Douglas and Gunn scheme appears to have an advantage over locally 
one-dimensional ( LOD ) or splitting schemes in which the intermediate steps do not 
satisfy the consistency condition [ 2 ]. As pointed out by Douglas and Gunn [ 12 ] the 
mixed derivatives can also be treated implicitly within the ADI framework, but this would 
increase the number of intermediate steps and thereby complicate the solution procedure. 

The above alternating-directional implicit technique, is applied to the conserva- 
tive equations. Then the equations are transformed into a now set of coordinates ( ?/, Q ) 
as discussed in section 2.8. Discretizing the spatial derivatives according to the three point 
central difference approximations for both second and first order derivatives, the following 
equations are obtained. 


Continuity equation : 
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x-momentum equation : 
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Energy equation : 
After first ADI step 
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3-4 Discretised Boundary Conditions 

Along the wall, C = 1 ; 
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The boundary condition for density is obtained as follows. 
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After expanding the nonlinear terms, transforming and discretizing the above equation 

becomes 
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Along the downstream boundary : y = 
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For density 
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It can be observed from the above set of discretized equations that, during first 
ADI step, the. continuity equation, ai-momentum equation and energy equations are only 
coupled and the ij moi 0.entum equation is decoupled from the other equations. Similarly, 
duiing second ADI step the continuity equation, the //-momentum equation and energy 

equations aie coupled and the ^-momentum equation is decoupled from the other equa- 
tions. 


Dming the first ADI step, the coupled continuity equation , x-momentiun equa- 
tion and energy equations along with corresponding boundary conditions for p,'u, and T 
foim a 3 x 3 block tridiagonal algebraic equations and can be solved for \E r * , T* and 
[ ^ ]• ^ ie decoupled y-momentum equation and the boundary condition for v form a 
simple tridiagonal algebraic equations and can be solved for '£* . 

Similarly, during the second ADI step, , ip" 11 and 'P^ +1 can be obtained 

by solving a set of coupi e( j equations, and \&“ +1 can be obtained independently. 

Fuithci, the solution at n + 1 time level is used for updating the turbulent 
viscosity. This piocedur e i s repeated until the steady state solution is achieved. 



Chapter 4 


Results and Discussion 


The numerical scheme as described in the earlier chapters has been implemented to analyse 
the flow over the tliin flat plate by solving the ensemble-averaged Navier-Stokes equations. 
A computer code has been developed and tested with the available experimental data for 
both laminar and turbulent regions. 

Firstly the laminar flow is analysed for a low temperature, low Mach no. subsonic 
flow. Such a low Mach no. flow behaves almost, incompressible and the code prediction 
can be compared with the Blasius solution. A plot of u vs ?/ & ( where r/ is equal to 
as in Blasius solution) is shown in Figure 4.1. The v vs 77 and T vs ?/ & are plotted in Figure 
4.2 and Figure 4.3 respectively. The coefficient of skin friction Cj vs Ra x is plotted in 
Figure 4.4. There is little variation in temperature T and velocity v and can be observed 
from the corresponding plots. Some minor fluctuations within 2 to 3 % are observed in 
the results and can be attributed to numerical oscillations. The variation of density is very 
little and can be negligible. The matching of the curves a vs 1 / and Cr vs Rc x with the 
known Blasius solution validate the code. The input parameters for this comparison of 
laminar flow are as follows. 


L 

= 1.0; 

'Umax 

= 0.125 

N x 

= 21 

Ny 

= 21 

Vx 

= 0.99 ; 

a y 

= 0.001 

M 

= 0.0718 ; 

Ret 

= 3.3 X 10 5 

Pr 

= 0.72 




The value of a x is taken as 0.99, which produces almost uniform grids in x- 
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direction, whereas a y is taken as 0.001, which cluster the grid points near the wall in 
.(/-direction as shown in Figure 2.2. The time step is chosen initially such that the CFL 
no. is equal to 2 and gradually changed to 20. Further higher CFL nos. lead to numerical 
oscillations in the solutions. The computation is stopped by checking the difference between 
two consecutive solutions at grid point (Nx/2,4 ) till it becomes less than the stipulated 
accuracy. The total number of time steps to get the steady state within the accuracy of 
0.000001 is found to be 180 for the above set of input parameters. 

There can be some error in the solution near the upstream boundary and near 
the down stream boundary, due to the leading edge singularity and to the' modelling the 
down stream boundary conditions. About 4 grid points on the either end, the result can 
be suspected to be erroneous to some extent. 

While analysing the laminar flow field, the turbulent viscosity is set to zero. It 
is observed that there is no need of the artificial viscosity and fourth order damping for 
the laminar flow. 

Next the turbulent region is analysed starting from leading edge ( including the 
laminar and transition regions ). A zero-equation turbulence model developed by Cebeci 
and Smith [ 8 ] is used for the turbulent viscosity in the turbulent zone. The transition is 
considered by taking an empirical expression for intermittency factor as discussed earlier in 
section 2.6. The code prediction is compared with the Simpson’s [ 24 ] experimental data 
for the mean velocity profiles with the following set of input parameters. 


L 

= 1.0; 

Umax 

= 0.055 

N x 

= 51 

Ny 

= 51 

<?x 

= 0.99; 


= 0.001 

M 

= 0.038 ; 

Re L 

= 2.0 x 10' 

Pr 

= 0.72 




The transition is assumed to occur at Re x = 4 x 10 5 [ 24 ], the CFL no. 

of about 10 is achieved for this calculation starting with 2 initially. The chock for the 
steady state is made as explained as earlier. Th('. result along with the Simpson’s [ 24 ] 
experimental data are shown in the Figure 4.5. This comparison validates the code for the 
turbulent flow region. 

The code is then used for 3 sets of input parameters to see the behaviour of code 
prediction. The various input parameters for these test, runs are given below. 
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Parameter 

Run 1 

Run 2 

Run 3 

L 

1 

1 

1 

Umax 

0.055 

0.055 

0.055 

Grid points 

51 x 51 

51 x 51 

101 x 101 

Rc l 

7.5 x 10 6 

6.0 x 10 (i 

1.0 x 10 7 

M 

0.36 

0.28 

0.5 

Pr 

0.72 

0.72 

0.72 


The Re xtr no. is taken as 4 x 10 6 as before for all the cases. To develop the 
solution, the CFL no. is initially taken as 2 and then gradually increased to 10. The total 
no. of time steps are found to be about 600 to get steady state. 

The Cf vs Re x for ail these 3 runs are plotted in Figure 4.6 along with the 
empirical correlation given in Schlichting [ 22 ] and log Cj vs log Re x plot is given in 
Figure 4.7. The correlation is restricted to a certain range of Reynolds no. but without any 
specification of the transition Reynolds no, Mach no., and Prandtl no.. However in spite of 
some minor differences between the predicted and empirical relation of the Cj with Re.,., in 
view of the earlier comparison of velocity profile with Simpson’s [ 24 ] experimental data, 
it can be concluded that the code predictions are reasonably accurate. 

The variation of mean velocity u with y for various Re x numbers ( at various 
locations on the plate for a single run, Run 3 ) are shown in Figure 4.8. Similarly the 
variations in mean temperature T and mean normal velocity v with y are plotted in Figures 
4.9. and 4.10 respectively. There occurs rise in temperature in the downstream direction 
due to viscous effects and can be seen in the Figure 4.9. The mean velocity v is becoming 
zero as the upstream boundary approaches. The minor oscillations can be attributed to 
numerical errors. 

The behaviour of the turbulence model, ie. the variation of the turbulent vis- 
cosity for Run 3, at Re x = 5 x 10 6 is plotted in Figure. 4 .11. Until the point A the 

turbulent viscosity is obtained from the inner layer expression and then from the A to B 
the outer layer expression is used as described in section 2.6. The convergence history of 
the velocity u with time at a grid point (51,4) for the Run 3 is shown in the Figure 4.12. In 
order to develop the mean velocity gradients, the turbulent viscosity model is introduced 
after some time steps ( presently, after 50 time steps ). A sudden rise in the velocity can 
be observed due to this. 




Figure 4.1: Comparison with Blasius Solution. 



Figure 4.2: Temperature Profiles for Laminar Case. 
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Figure 4.4: Local Skin Friction Coefficient vs Re x f° r Laminar case 
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Figure 4.6: Local Skin Friction Coefficient, vs Re x For Turbulent Flow 
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Figure 4.8: Mean u Velocity Profiles at Different Re, r For Turbulent Flow 
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4.9 Mean T Temperature Profiles at Different Re- For turbulent flow. 
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4.10 Mean v Velocity Profiles at Different Re,. For Turbulent Flow. 
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Figure 4.11: Turbulent Viscosity vs y at Re 
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Figure 4.12: Convergence History of u at Grid Point (51,4) for Run 3 






Chapt er 5 


Conclusions and Suggested 
Extensions 

5.1 Conclusions 

The time dependent, multi-dimensional, implicit finite difference scheme has been suc- 
cessfully applied to solve ensemble-averaged Navier- Stokes equations for a compressible, 
turbulent flow over a thin flat plate . The validity of the associated computer code has 
been established by comparing the code prediction with experimental/calculated results. 
A second order artificial viscosity model is used for the high Reynolds no. flows without 
much loss in accuracy in the solution to suppress the numerical oscillations. A fourth order 
dissipation term is used to smoothen the solution. 

In the laminar region, the calculated velocity profile agreed with the well-known 
Blasius solution. For the turbulent flow, zero-equation Cebeci and Smith model is used 
and a reasonable agreement is found between the code predicted mean velocity u profile 
with the Simpson’s experimental data. The calculated coefficient of local skin friction is 
compared with the known empirical correlation for it. Three computer runs are made for 
different input parameters and the results are plotted. 



50 


5.2 Suggested Extensions 


The present, numerical scheme and the associated computer code is capable of handling 
two-dimensional , non-reacting turbulent flow fields. It can be extended to deal with three- 
dimensional flows such as in rectangular intakes etc. including higher order turbulence 

models. 


Another extension can be made with multicomponent , reacting flows winch could 
be used, for instance, to solve the problem of the fiat solid propellent grain burninig 
( Emmon’s flat plate burning problem ). 
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